module UrbanFluxesMod

  !----------------------------------------------------------------------- 
  ! !DESCRIPTION: 
  ! Calculate solar and longwave radiation, and turbulent fluxes for urban landunit
  !
  ! !USES:
  use shr_kind_mod         , only : r8 => shr_kind_r8
  use shr_sys_mod          , only : shr_sys_flush 
  use shr_log_mod          , only : errMsg => shr_log_errMsg
  use decompMod            , only : bounds_type, subgrid_level_landunit
  use clm_varpar           , only : numrad
  use clm_varctl           , only : iulog
  use abortutils           , only : endrun  
  use UrbanParamsType      , only : urbanparams_type
  use UrbanParamsType      , only : urban_wasteheat_on, urban_hac_on, urban_hac
  use UrbanParamsType      , only : IsSimpleBuildTemp
  use atm2lndType          , only : atm2lnd_type
  use SoilStateType        , only : soilstate_type
  use TemperatureType      , only : temperature_type
  use WaterStateBulkType       , only : waterstatebulk_type
  use WaterDiagnosticBulkType       , only : waterdiagnosticbulk_type
  use FrictionVelocityMod  , only : frictionvel_type
  use EnergyFluxType       , only : energyflux_type
  use WaterFluxBulkType        , only : waterfluxbulk_type
  use Wateratm2lndBulkType        , only : wateratm2lndbulk_type
  use HumanIndexMod        , only : humanindex_type
  use GridcellType         , only : grc                
  use LandunitType         , only : lun                
  use ColumnType           , only : col                
  use PatchType            , only : patch                
  !
  ! !PUBLIC TYPES:
  implicit none
  save
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  public :: UrbanFluxes       ! Urban physics - turbulent fluxes
  public :: readParams
  !-----------------------------------------------------------------------

  ! !PRIVATE FUNCTIONS:
  private :: wasteheat               ! Figure out the energy flux from urban heating and cooling
  private :: simple_wasteheatfromac  ! Calculate waste heat from air-conditioning with the simpler method (CLM4.5)
  private :: calc_simple_internal_building_temp ! Calculate internal building temperature by simpler method (CLM4.5)

  type, private :: params_type
     real(r8) :: wind_min ! Minimum wind speed at the atmospheric forcing height (m/s)
  end type params_type
  type(params_type), private ::  params_inst

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

contains

  !------------------------------------------------------------------------------
  subroutine readParams( ncid )
    !
    ! !USES:
    use ncdio_pio, only: file_desc_t
    use paramUtilMod, only: readNcdioScalar
    !
    ! !ARGUMENTS:
    implicit none
    type(file_desc_t),intent(inout) :: ncid   ! pio netCDF file id
    !
    ! !LOCAL VARIABLES:
    character(len=*), parameter :: subname = 'readParams_UrbanFluxes'
    !--------------------------------------------------------------------

    ! Minimum wind speed at the atmospheric forcing height (m/s)
    call readNcdioScalar(ncid, 'wind_min', subname, params_inst%wind_min)

  end subroutine readParams

  !-----------------------------------------------------------------------
  subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl,                        &
       num_urbanl, filter_urbanl, num_urbanc, filter_urbanc, num_urbanp, filter_urbanp, &
       atm2lnd_inst, urbanparams_inst, soilstate_inst, temperature_inst,                &
       waterstatebulk_inst, waterdiagnosticbulk_inst, frictionvel_inst,                 &
       energyflux_inst, waterfluxbulk_inst, wateratm2lndbulk_inst,                      &
       humanindex_inst) 
    !
    ! !DESCRIPTION: 
    ! Turbulent and momentum fluxes from urban canyon (consisting of roof, sunwall, 
    ! shadewall, pervious and impervious road).

    ! !USES:
    use clm_varcon          , only : cpair, vkc, spval, grav, pondmx_urban, rpi, rgas
    use clm_varcon          , only : ht_wasteheat_factor, ac_wasteheat_factor, wasteheat_limit
    use column_varcon       , only : icol_shadewall, icol_road_perv, icol_road_imperv
    use column_varcon       , only : icol_roof, icol_sunwall
    use filterMod           , only : filter
    use QSatMod             , only : QSat
    use clm_varpar          , only : maxpatch_urb, nlevurb, nlevgrnd
    use clm_time_manager    , only : get_step_size_real, get_nstep, is_near_local_noon
    use HumanIndexMod       , only : all_human_stress_indices, fast_human_stress_indices, &
                                     Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, &
                                     swbgt, hmdex, dis_coi, dis_coiS, THIndex, &
                                     SwampCoolEff, KtoC, VaporPres
    !
    ! !ARGUMENTS:
    type(bounds_type)      , intent(in)    :: bounds    
    integer                , intent(in)    :: num_nourbanl       ! number of non-urban landunits in clump
    integer                , intent(in)    :: filter_nourbanl(:) ! non-urban landunit filter
    integer                , intent(in)    :: num_urbanl         ! number of urban landunits in clump
    integer                , intent(in)    :: filter_urbanl(:)   ! urban landunit filter
    integer                , intent(in)    :: num_urbanc         ! number of urban columns in clump
    integer                , intent(in)    :: filter_urbanc(:)   ! urban column filter
    integer                , intent(in)    :: num_urbanp         ! number of urban patches in clump
    integer                , intent(in)    :: filter_urbanp(:)   ! urban pft filter
    type(atm2lnd_type)     , intent(in)    :: atm2lnd_inst
    type(urbanparams_type) , intent(in)    :: urbanparams_inst
    type(soilstate_type)   , intent(inout) :: soilstate_inst
    type(temperature_type) , intent(inout) :: temperature_inst
    type(waterstatebulk_type)  , intent(inout) :: waterstatebulk_inst
    type(waterdiagnosticbulk_type)  , intent(inout) :: waterdiagnosticbulk_inst
    type(frictionvel_type) , intent(inout) :: frictionvel_inst
    type(waterfluxbulk_type)   , intent(inout) :: waterfluxbulk_inst
    type(wateratm2lndbulk_type)   , intent(inout) :: wateratm2lndbulk_inst
    type(energyflux_type)  , intent(inout) :: energyflux_inst
    type(humanindex_type)  , intent(inout) :: humanindex_inst
    !
    ! !LOCAL VARIABLES:
    character(len=*), parameter :: sub="UrbanFluxes"
    integer  :: fp,fc,fl,f,p,c,l,g,j,pi,i     ! indices

    real(r8) :: canyontop_wind(bounds%begl:bounds%endl)              ! wind at canyon top (m/s) 
    real(r8) :: canyon_u_wind(bounds%begl:bounds%endl)               ! u-component of wind speed inside canyon (m/s)
    real(r8) :: canyon_wind(bounds%begl:bounds%endl)                 ! net wind speed inside canyon (m/s)
    real(r8) :: canyon_resistance(bounds%begl:bounds%endl)           ! resistance to heat and moisture transfer from canyon road/walls to canyon air (s/m)

    real(r8) :: ur(bounds%begl:bounds%endl)                          ! wind speed at reference height (m/s)
    real(r8) :: ustar(bounds%begl:bounds%endl)                       ! friction velocity (m/s)
    real(r8) :: ramu(bounds%begl:bounds%endl)                        ! aerodynamic resistance (s/m)
    real(r8) :: rahu(bounds%begl:bounds%endl)                        ! thermal resistance (s/m)
    real(r8) :: rawu(bounds%begl:bounds%endl)                        ! moisture resistance (s/m)
    real(r8) :: temp1(bounds%begl:bounds%endl)                       ! relation for potential temperature profile
    real(r8) :: temp12m(bounds%begl:bounds%endl)                     ! relation for potential temperature profile applied at 2-m
    real(r8) :: temp2(bounds%begl:bounds%endl)                       ! relation for specific humidity profile
    real(r8) :: temp22m(bounds%begl:bounds%endl)                     ! relation for specific humidity profile applied at 2-m
    real(r8) :: thm_g(bounds%begl:bounds%endl)                       ! intermediate variable (forc_t+0.0098*forc_hgt_t)
    real(r8) :: thv_g(bounds%begl:bounds%endl)                       ! virtual potential temperature (K)
    real(r8) :: dth(bounds%begl:bounds%endl)                         ! diff of virtual temp. between ref. height and surface
    real(r8) :: dqh(bounds%begl:bounds%endl)                         ! diff of humidity between ref. height and surface
    real(r8) :: zldis(bounds%begl:bounds%endl)                       ! reference height "minus" zero displacement height (m)
    real(r8) :: um(bounds%begl:bounds%endl)                          ! wind speed including the stablity effect (m/s)
    real(r8) :: obu(bounds%begl:bounds%endl)                         ! Monin-Obukhov length (m)
    real(r8) :: taf_numer(bounds%begl:bounds%endl)                   ! numerator of taf equation (K m/s)
    real(r8) :: taf_denom(bounds%begl:bounds%endl)                   ! denominator of taf equation (m/s)
    real(r8) :: qaf_numer(bounds%begl:bounds%endl)                   ! numerator of qaf equation (kg m/kg s)
    real(r8) :: qaf_denom(bounds%begl:bounds%endl)                   ! denominator of qaf equation (m/s)
    real(r8) :: wtas(bounds%begl:bounds%endl)                        ! sensible heat conductance for urban air to atmospheric air (m/s)
    real(r8) :: wtaq(bounds%begl:bounds%endl)                        ! latent heat conductance for urban air to atmospheric air (m/s)
    real(r8) :: wts_sum(bounds%begl:bounds%endl)                     ! sum of wtas, wtus_roof, wtus_road_perv, wtus_road_imperv, wtus_sunwall, wtus_shadewall
    real(r8) :: wtq_sum(bounds%begl:bounds%endl)                     ! sum of wtaq, wtuq_roof, wtuq_road_perv, wtuq_road_imperv, wtuq_sunwall, wtuq_shadewall
    real(r8) :: beta(bounds%begl:bounds%endl)                        ! coefficient of convective velocity
    real(r8) :: zii(bounds%begl:bounds%endl)                         ! convective boundary layer height (m)
    real(r8) :: fm(bounds%begl:bounds%endl)                          ! needed for BGC only to diagnose 10m wind speed
    real(r8) :: wtus(bounds%begc:bounds%endc)                        ! sensible heat conductance for urban columns (scaled) (m/s)
    real(r8) :: wtuq(bounds%begc:bounds%endc)                        ! latent heat conductance for urban columns (scaled) (m/s)
    integer  :: iter                                                 ! iteration index
    real(r8) :: dthv                                                 ! diff of vir. poten. temp. between ref. height and surface
    real(r8) :: tstar                                                ! temperature scaling parameter
    real(r8) :: qstar                                                ! moisture scaling parameter
    real(r8) :: thvstar                                              ! virtual potential temperature scaling parameter
    real(r8) :: wtus_roof(bounds%begl:bounds%endl)                   ! sensible heat conductance for roof (scaled) (m/s)
    real(r8) :: wtuq_roof(bounds%begl:bounds%endl)                   ! latent heat conductance for roof (scaled) (m/s)
    real(r8) :: wtus_road_perv(bounds%begl:bounds%endl)              ! sensible heat conductance for pervious road (scaled) (m/s)
    real(r8) :: wtuq_road_perv(bounds%begl:bounds%endl)              ! latent heat conductance for pervious road (scaled) (m/s)
    real(r8) :: wtus_road_imperv(bounds%begl:bounds%endl)            ! sensible heat conductance for impervious road (scaled) (m/s)
    real(r8) :: wtuq_road_imperv(bounds%begl:bounds%endl)            ! latent heat conductance for impervious road (scaled) (m/s)
    real(r8) :: wtus_sunwall(bounds%begl:bounds%endl)                ! sensible heat conductance for sunwall (scaled) (m/s)
    real(r8) :: wtuq_sunwall(bounds%begl:bounds%endl)                ! latent heat conductance for sunwall (scaled) (m/s)
    real(r8) :: wtus_shadewall(bounds%begl:bounds%endl)              ! sensible heat conductance for shadewall (scaled) (m/s)
    real(r8) :: wtuq_shadewall(bounds%begl:bounds%endl)              ! latent heat conductance for shadewall (scaled) (m/s)
    real(r8) :: wtus_roof_unscl(bounds%begl:bounds%endl)             ! sensible heat conductance for roof (not scaled) (m/s)
    real(r8) :: wtuq_roof_unscl(bounds%begl:bounds%endl)             ! latent heat conductance for roof (not scaled) (m/s)
    real(r8) :: wtus_road_perv_unscl(bounds%begl:bounds%endl)        ! sensible heat conductance for pervious road (not scaled) (m/s)
    real(r8) :: wtuq_road_perv_unscl(bounds%begl:bounds%endl)        ! latent heat conductance for pervious road (not scaled) (m/s)
    real(r8) :: wtus_road_imperv_unscl(bounds%begl:bounds%endl)      ! sensible heat conductance for impervious road (not scaled) (m/s)
    real(r8) :: wtuq_road_imperv_unscl(bounds%begl:bounds%endl)      ! latent heat conductance for impervious road (not scaled) (m/s)
    real(r8) :: wtus_sunwall_unscl(bounds%begl:bounds%endl)          ! sensible heat conductance for sunwall (not scaled) (m/s)
    real(r8) :: wtuq_sunwall_unscl(bounds%begl:bounds%endl)          ! latent heat conductance for sunwall (not scaled) (m/s)
    real(r8) :: wtus_shadewall_unscl(bounds%begl:bounds%endl)        ! sensible heat conductance for shadewall (not scaled) (m/s)
    real(r8) :: wtuq_shadewall_unscl(bounds%begl:bounds%endl)        ! latent heat conductance for shadewall (not scaled) (m/s)
    real(r8) :: wc                                                   ! convective velocity (m/s)
    real(r8) :: zeta                                                 ! dimensionless height used in Monin-Obukhov theory 
    real(r8) :: eflx_sh_grnd_scale(bounds%begp:bounds%endp)          ! scaled sensible heat flux from ground (W/m**2) [+ to atm] 
    real(r8) :: qflx_evap_soi_scale(bounds%begp:bounds%endp)         ! scaled soil evaporation (mm H2O/s) (+ = to atm) 
    real(r8) :: eflx_wasteheat_roof(bounds%begl:bounds%endl)         ! sensible heat flux from urban heating/cooling sources of waste heat for roof (W/m**2)
    real(r8) :: eflx_wasteheat_sunwall(bounds%begl:bounds%endl)      ! sensible heat flux from urban heating/cooling sources of waste heat for sunwall (W/m**2)
    real(r8) :: eflx_wasteheat_shadewall(bounds%begl:bounds%endl)    ! sensible heat flux from urban heating/cooling sources of waste heat for shadewall (W/m**2)
    real(r8) :: eflx_heat_from_ac_roof(bounds%begl:bounds%endl)      ! sensible heat flux put back into canyon due to heat removal by AC for roof (W/m**2)
    real(r8) :: eflx_heat_from_ac_sunwall(bounds%begl:bounds%endl)   ! sensible heat flux put back into canyon due to heat removal by AC for sunwall (W/m**2)
    real(r8) :: eflx_heat_from_ac_shadewall(bounds%begl:bounds%endl) ! sensible heat flux put back into canyon due to heat removal by AC for shadewall (W/m**2)
    real(r8) :: eflx(bounds%begl:bounds%endl)                        ! total sensible heat flux for error check (W/m**2)
    real(r8) :: qflx(bounds%begl:bounds%endl)                        ! total water vapor flux for error check (kg/m**2/s)
    real(r8) :: eflx_scale(bounds%begl:bounds%endl)                  ! sum of scaled sensible heat fluxes for urban columns for error check (W/m**2)
    real(r8) :: qflx_scale(bounds%begl:bounds%endl)                  ! sum of scaled water vapor fluxes for urban columns for error check (kg/m**2/s)
    real(r8) :: eflx_err(bounds%begl:bounds%endl)                    ! sensible heat flux error (W/m**2)
    real(r8) :: qflx_err(bounds%begl:bounds%endl)                    ! water vapor flux error (kg/m**2/s)
    real(r8) :: fwet_roof                                            ! fraction of roof surface that is wet (-)
    real(r8) :: fwet_road_imperv                                     ! fraction of impervious road surface that is wet (-)
    real(r8) :: dtime                                                ! land model time step (sec)
    logical  :: found                                                ! flag in search loop
    integer  :: indexl                                               ! index of first found in search loop
    integer  :: nstep                                                ! time step number
    real(r8) :: e_ref2m                                              ! 2 m height surface saturated vapor pressure [Pa]
    real(r8) :: qsat_ref2m                                           ! 2 m height surface saturated specific humidity [kg/kg]
    !
    real(r8), parameter :: lapse_rate = 0.0098_r8 ! Dry adiabatic lapse rate (K/m)
    integer , parameter  :: niters = 3            ! maximum number of iterations for surface temperature
    !-----------------------------------------------------------------------

    associate(                                                                & 
         snl                 =>   col%snl                                   , & ! Input:  [integer  (:)   ]  number of snow layers                              
         ctype               =>   col%itype                                 , & ! Input:  [integer  (:)   ]  column type                                        
         z_0_town            =>   lun%z_0_town                              , & ! Input:  [real(r8) (:)   ]  momentum roughness length of urban landunit (m)   
         z_d_town            =>   lun%z_d_town                              , & ! Input:  [real(r8) (:)   ]  displacement height of urban landunit (m)         
         ht_roof             =>   lun%ht_roof                               , & ! Input:  [real(r8) (:)   ]  height of urban roof (m)                          
         wtlunit_roof        =>   lun%wtlunit_roof                          , & ! Input:  [real(r8) (:)   ]  weight of roof with respect to landunit           
         canyon_hwr          =>   lun%canyon_hwr                            , & ! Input:  [real(r8) (:)   ]  ratio of building height to street width          
         wtroad_perv         =>   lun%wtroad_perv                           , & ! Input:  [real(r8) (:)   ]  weight of pervious road wrt total road            

         forc_t              =>   atm2lnd_inst%forc_t_not_downscaled_grc    , & ! Input:  [real(r8) (:)   ]  atmospheric temperature (K)                       
         forc_th             =>   atm2lnd_inst%forc_th_not_downscaled_grc   , & ! Input:  [real(r8) (:)   ]  atmospheric potential temperature (K)             
         forc_rho            =>   atm2lnd_inst%forc_rho_not_downscaled_grc  , & ! Input:  [real(r8) (:)   ]  density (kg/m**3)                                 
         forc_q              =>   wateratm2lndbulk_inst%forc_q_not_downscaled_grc    , & ! Input:  [real(r8) (:)   ]  atmospheric specific humidity (kg/kg)             
         forc_pbot           =>   atm2lnd_inst%forc_pbot_not_downscaled_grc , & ! Input:  [real(r8) (:)   ]  atmospheric pressure (Pa)                         
         forc_u              =>   atm2lnd_inst%forc_u_grc                   , & ! Input:  [real(r8) (:)   ]  atmospheric wind speed in east direction (m/s)    
         forc_v              =>   atm2lnd_inst%forc_v_grc                   , & ! Input:  [real(r8) (:)   ]  atmospheric wind speed in north direction (m/s)   

         wind_hgt_canyon     =>   urbanparams_inst%wind_hgt_canyon          , & ! Input:  [real(r8) (:)   ]  height above road at which wind in canyon is to be computed (m)
         eflx_traffic_factor =>   urbanparams_inst%eflx_traffic_factor      , & ! Input:  [real(r8) (:)   ]  multiplicative urban traffic factor for sensible heat flux

         rootr_road_perv     =>   soilstate_inst%rootr_road_perv_col        , & ! Input:  [real(r8) (:,:) ]  effective fraction of roots in each soil layer for urban pervious road
         soilalpha_u         =>   soilstate_inst%soilalpha_u_col            , & ! Input:  [real(r8) (:)   ]  Urban factor that reduces ground saturated specific humidity (-)
         rootr               =>   soilstate_inst%rootr_patch                , & ! Output: [real(r8) (:,:) ]  effective fraction of roots in each soil layer (SMS method only) 

         t_grnd              =>   temperature_inst%t_grnd_col               , & ! Input:  [real(r8) (:)   ]  ground surface temperature (K)                    
         t_soisno            =>   temperature_inst%t_soisno_col             , & ! Input:  [real(r8) (:,:) ]  soil temperature (K)                            
         t_ref2m             =>   temperature_inst%t_ref2m_patch            , & ! Output: [real(r8) (:)   ]  2 m height surface air temperature (K)            
         t_ref2m_u           =>   temperature_inst%t_ref2m_u_patch          , & ! Output: [real(r8) (:)   ]  Urban 2 m height surface air temperature (K)     
         t_veg               =>   temperature_inst%t_veg_patch              , & ! Output: [real(r8) (:)   ]  vegetation temperature (K)                        
         taf                 =>   temperature_inst%taf_lun                  , & ! Output: [real(r8) (:)   ]  urban canopy air temperature (K)                  


         tc_ref2m            => humanindex_inst%tc_ref2m_patch              , & ! Output: [real(r8) (:)   ]  2 m height surface air temperature (C)
         vap_ref2m           => humanindex_inst%vap_ref2m_patch             , & ! Output: [real(r8) (:)   ]  2 m height vapor pressure (Pa)
         appar_temp_ref2m    => humanindex_inst%appar_temp_ref2m_patch      , & ! Output: [real(r8) (:)   ]  2 m apparent temperature (C)
         appar_temp_ref2m_u  => humanindex_inst%appar_temp_ref2m_u_patch    , & ! Output: [real(r8) (:)   ]  Urban 2 m apparent temperature (C)
         swbgt_ref2m         => humanindex_inst%swbgt_ref2m_patch           , & ! Output: [real(r8) (:)   ]  2 m Simplified Wetbulb Globe temperature (C)
         swbgt_ref2m_u       => humanindex_inst%swbgt_ref2m_u_patch         , & ! Output: [real(r8) (:)   ]  Urban 2 m Simplified Wetbulb Globe temperature (C)
         humidex_ref2m       => humanindex_inst%humidex_ref2m_patch         , & ! Output: [real(r8) (:)   ]  2 m Humidex (C)
         humidex_ref2m_u     => humanindex_inst%humidex_ref2m_u_patch       , & ! Output: [real(r8) (:)   ]  Urban 2 m Humidex (C)
         wbt_ref2m           => humanindex_inst%wbt_ref2m_patch             , & ! Output: [real(r8) (:)   ]  2 m Stull Wet Bulb temperature (C)
         wbt_ref2m_u         => humanindex_inst%wbt_ref2m_u_patch           , & ! Output: [real(r8) (:)   ]  Urban 2 m Stull Wet Bulb temperature (C)
         wb_ref2m            => humanindex_inst%wb_ref2m_patch              , & ! Output: [real(r8) (:)   ]  2 m Wet Bulb temperature (C)
         wb_ref2m_u          => humanindex_inst%wb_ref2m_u_patch            , & ! Output: [real(r8) (:)   ]  Urban 2 m Wet Bulb temperature (C)
         teq_ref2m           => humanindex_inst%teq_ref2m_patch             , & ! Output: [real(r8) (:)   ]  2 m height Equivalent temperature (K)
         teq_ref2m_u         => humanindex_inst%teq_ref2m_u_patch           , & ! Output: [real(r8) (:)   ]  Urban 2 m Equivalent temperature (K)
         ept_ref2m           => humanindex_inst%ept_ref2m_patch             , & ! Output: [real(r8) (:)   ]  2 m height Equivalent Potential temperature (K)
         ept_ref2m_u         => humanindex_inst%ept_ref2m_u_patch           , & ! Output: [real(r8) (:)   ]  Urban 2 m height Equivalent Potential temperature (K)
         discomf_index_ref2m     => humanindex_inst%discomf_index_ref2m_patch   , & ! Output: [real(r8) (:)   ]  2 m Discomfort Index temperature (C)
         discomf_index_ref2m_u   => humanindex_inst%discomf_index_ref2m_u_patch , & ! Output: [real(r8) (:)   ]  Urban 2 m Discomfort Index temperature (C)
         discomf_index_ref2mS    => humanindex_inst%discomf_index_ref2mS_patch  , & ! Output: [real(r8) (:)   ]  2 m height Discomfort Index Stull temperature (C)
         discomf_index_ref2mS_u  => humanindex_inst%discomf_index_ref2mS_u_patch, & ! Output: [real(r8) (:)   ]  Urban 2 m Discomfort Index Stull temperature (K)
         nws_hi_ref2m        => humanindex_inst%nws_hi_ref2m_patch          , & ! Output: [real(r8) (:)   ]  2 m NWS Heat Index (C)
         nws_hi_ref2m_u      => humanindex_inst%nws_hi_ref2m_u_patch        , & ! Output: [real(r8) (:)   ]  Urban 2 m NWS Heat Index (C)
         thip_ref2m          => humanindex_inst%thip_ref2m_patch            , & ! Output: [real(r8) (:)   ]  2 m Temperature Humidity Index Physiology (C)
         thip_ref2m_u        => humanindex_inst%thip_ref2m_u_patch          , & ! Output: [real(r8) (:)   ]  Urban 2 m Temperature Humidity Index Physiology (C)
         thic_ref2m          => humanindex_inst%thic_ref2m_patch            , & ! Output: [real(r8) (:)   ]  2 m Temperature Humidity Index Comfort (C)
         thic_ref2m_u        => humanindex_inst%thic_ref2m_u_patch          , & ! Output: [real(r8) (:)   ]  Urban 2 m Temperature Humidity Index Comfort (C)
         swmp65_ref2m        => humanindex_inst%swmp65_ref2m_patch          , & ! Output: [real(r8) (:)   ]  2 m Swamp Cooler temperature 65% effi (C)
         swmp65_ref2m_u      => humanindex_inst%swmp65_ref2m_u_patch        , & ! Output: [real(r8) (:)   ]  Urban 2 m Swamp Cooler temperature 65% effi (C)
         swmp80_ref2m        => humanindex_inst%swmp80_ref2m_patch          , & ! Output: [real(r8) (:)   ]  2 m Swamp Cooler temperature 80% effi (C)
         swmp80_ref2m_u      => humanindex_inst%swmp80_ref2m_u_patch        , & ! Output: [real(r8) (:)   ]  Urban 2 m Swamp Cooler temperature 80% effi (C)

         frac_sno            =>   waterdiagnosticbulk_inst%frac_sno_col              , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by snow (0 to 1)       
         snow_depth          =>   waterdiagnosticbulk_inst%snow_depth_col            , & ! Input:  [real(r8) (:)   ]  snow height (m)                                   
         dqgdT               =>   waterdiagnosticbulk_inst%dqgdT_col                 , & ! Input:  [real(r8) (:)   ]  temperature derivative of "qg"                    
         qg                  =>   waterdiagnosticbulk_inst%qg_col                    , & ! Input:  [real(r8) (:)   ]  specific humidity at ground surface (kg/kg)       
         h2osoi_ice          =>   waterstatebulk_inst%h2osoi_ice_col            , & ! Input:  [real(r8) (:,:) ]  ice lens (kg/m2)                                
         h2osoi_liq          =>   waterstatebulk_inst%h2osoi_liq_col            , & ! Input:  [real(r8) (:,:) ]  liquid water (kg/m2)                            
         qaf                 =>   waterdiagnosticbulk_inst%qaf_lun                   , & ! Output: [real(r8) (:)   ]  urban canopy air specific humidity (kg/kg)        
         q_ref2m             =>   waterdiagnosticbulk_inst%q_ref2m_patch             , & ! Output: [real(r8) (:)   ]  2 m height surface specific humidity (kg/kg)      
         rh_ref2m            =>   waterdiagnosticbulk_inst%rh_ref2m_patch            , & ! Output: [real(r8) (:)   ]  2 m height surface relative humidity (%)          
         rh_ref2m_u          =>   waterdiagnosticbulk_inst%rh_ref2m_u_patch          , & ! Output: [real(r8) (:)   ]  2 m height surface relative humidity (%)          

         forc_hgt_u_patch    =>   frictionvel_inst%forc_hgt_u_patch         , & ! Input:  [real(r8) (:)   ]  observational height of wind at patch-level (m)     
         forc_hgt_t_patch    =>   frictionvel_inst%forc_hgt_t_patch         , & ! Input:  [real(r8) (:)   ]  observational height of temperature at patch-level (m)
         zetamax             =>   frictionvel_inst%zetamaxstable            , & ! Input:  [real(r8)       ]  max zeta value under stable conditions
         ram1                =>   frictionvel_inst%ram1_patch               , & ! Output: [real(r8) (:)   ]  aerodynamical resistance (s/m)                    
         u10_clm             =>   frictionvel_inst%u10_clm_patch            , & ! Input:  [real(r8) (:)   ]  10 m height winds (m/s)

         htvp                =>   energyflux_inst%htvp_col                  , & ! Input:  [real(r8) (:)   ]  latent heat of evaporation (/sublimation) (J/kg)  
         dlrad               =>   energyflux_inst%dlrad_patch               , & ! Output: [real(r8) (:)   ]  downward longwave radiation below the canopy (W/m**2)
         ulrad               =>   energyflux_inst%ulrad_patch               , & ! Output: [real(r8) (:)   ]  upward longwave radiation above the canopy (W/m**2)
         cgrnds              =>   energyflux_inst%cgrnds_patch              , & ! Output: [real(r8) (:)   ]  deriv, of soil sensible heat flux wrt soil temp (W/m**2/K)
         cgrndl              =>   energyflux_inst%cgrndl_patch              , & ! Output: [real(r8) (:)   ]  deriv of soil latent heat flux wrt soil temp (W/m**2/K)
         cgrnd               =>   energyflux_inst%cgrnd_patch               , & ! Output: [real(r8) (:)   ]  deriv. of soil energy flux wrt to soil temp (W/m**2/K)
         eflx_sh_grnd        =>   energyflux_inst%eflx_sh_grnd_patch        , & ! Output: [real(r8) (:)   ]  sensible heat flux from ground (W/m**2) [+ to atm]
         eflx_sh_tot         =>   energyflux_inst%eflx_sh_tot_patch         , & ! Output: [real(r8) (:)   ]  total sensible heat flux (W/m**2) [+ to atm]      
         eflx_sh_tot_u       =>   energyflux_inst%eflx_sh_tot_u_patch       , & ! Output: [real(r8) (:)   ]  urban total sensible heat flux (W/m**2) [+ to atm]
         eflx_sh_snow        =>   energyflux_inst%eflx_sh_snow_patch        , & ! Output: [real(r8) (:)   ]  sensible heat flux from snow (W/m**2) [+ to atm]  
         eflx_sh_soil        =>   energyflux_inst%eflx_sh_soil_patch        , & ! Output: [real(r8) (:)   ]  sensible heat flux from soil (W/m**2) [+ to atm]  
         eflx_sh_h2osfc      =>   energyflux_inst%eflx_sh_h2osfc_patch      , & ! Output: [real(r8) (:)   ]  sensible heat flux from soil (W/m**2) [+ to atm]  
         eflx_traffic        =>   energyflux_inst%eflx_traffic_lun          , & ! Output: [real(r8) (:)   ]  traffic sensible heat flux (W/m**2)               
         eflx_wasteheat      =>   energyflux_inst%eflx_wasteheat_lun        , & ! Output: [real(r8) (:)   ]  sensible heat flux from urban heating/cooling sources of waste heat (W/m**2)
         eflx_urban_ac       =>   energyflux_inst%eflx_urban_ac_lun         , & ! Input:  [real(r8) (:)   ]  urban air conditioning flux (W/m**2)              
         eflx_heat_from_ac   =>   energyflux_inst%eflx_heat_from_ac_lun     , & ! Output: [real(r8) (:)   ]  sensible heat flux put back into canyon due to removal by AC (W/m**2)
         eflx_urban_heat     =>   energyflux_inst%eflx_urban_heat_lun       , & ! Input:  [real(r8) (:)   ]  urban heating flux (W/m**2)
         eflx_urban_ac_col   =>   energyflux_inst%eflx_urban_ac_col         , & ! Input:  [real(r8) (:)   ]  urban air conditioning flux (W/m**2)
         eflx_urban_heat_col =>   energyflux_inst%eflx_urban_heat_col       , & ! Input:  [real(r8) (:)   ]  urban heating flux (W/m**2)
         taux                =>   energyflux_inst%taux_patch                , & ! Output: [real(r8) (:)   ]  wind (shear) stress: e-w (kg/m/s**2)              
         tauy                =>   energyflux_inst%tauy_patch                , & ! Output: [real(r8) (:)   ]  wind (shear) stress: n-s (kg/m/s**2)               

         qflx_evap_soi       =>   waterfluxbulk_inst%qflx_evap_soi_patch        , & ! Output: [real(r8) (:)   ]  soil evaporation (mm H2O/s) (+ = to atm)          
         qflx_tran_veg       =>   waterfluxbulk_inst%qflx_tran_veg_patch        , & ! Output: [real(r8) (:)   ]  vegetation transpiration (mm H2O/s) (+ = to atm)  
         qflx_evap_veg       =>   waterfluxbulk_inst%qflx_evap_veg_patch        , & ! Output: [real(r8) (:)   ]  vegetation evaporation (mm H2O/s) (+ = to atm)    

         begl                =>   bounds%begl                               , &
         endl                =>   bounds%endl                                 &
         )

      ! Define fields that appear on the restart file for non-urban landunits 
      
      do fl = 1,num_nourbanl
         l = filter_nourbanl(fl)
         taf(l) = spval
         qaf(l) = spval
      end do

      ! Get time step
      nstep = get_nstep()

      ! Set constants (same as in Biogeophysics1Mod)
      beta(begl:endl) = 1._r8             ! Should be set to the same values as in Biogeophysics1Mod
      zii(begl:endl)  = 1000._r8          ! Should be set to the same values as in Biogeophysics1Mod

      ! Get current date
      dtime = get_step_size_real()

      ! Compute canyontop wind using Masson (2000)

      do fl = 1, num_urbanl
         l = filter_urbanl(fl)
         g = lun%gridcell(l)

         ! Error checks

         if (ht_roof(l) - z_d_town(l) <= z_0_town(l)) then
            write (iulog,*) 'aerodynamic parameter error in UrbanFluxes'
            write (iulog,*) 'h_r - z_d <= z_0'
            write (iulog,*) 'ht_roof, z_d_town, z_0_town: ', ht_roof(l), z_d_town(l), &
                 z_0_town(l)
            write (iulog,*) 'clm model is stopping'
            call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__))
         end if
         if (forc_hgt_u_patch(lun%patchi(l)) - z_d_town(l) <= z_0_town(l)) then
            write (iulog,*) 'aerodynamic parameter error in UrbanFluxes'
            write (iulog,*) 'h_u - z_d <= z_0'
            write (iulog,*) 'forc_hgt_u_patch, z_d_town, z_0_town: ', forc_hgt_u_patch(lun%patchi(l)), z_d_town(l), &
                 z_0_town(l)
            write (iulog,*) 'clm model is stopping'
            call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__))
         end if

         ! Magnitude of atmospheric wind

         ur(l) = max(params_inst%wind_min,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g)))

         ! Canyon top wind

         canyontop_wind(l) = ur(l) * &
              log( (ht_roof(l)-z_d_town(l)) / z_0_town(l) ) / &
              log( (forc_hgt_u_patch(lun%patchi(l))-z_d_town(l)) / z_0_town(l) )

         ! U component of canyon wind 

         if (canyon_hwr(l) < 0.5_r8) then  ! isolated roughness flow
            canyon_u_wind(l) = canyontop_wind(l) * exp( -0.5_r8*canyon_hwr(l)* &
                 (1._r8-(wind_hgt_canyon(l)/ht_roof(l))) )
         else if (canyon_hwr(l) < 1.0_r8) then ! wake interference flow
            canyon_u_wind(l) = canyontop_wind(l) * (1._r8+2._r8*(2._r8/rpi - 1._r8)* &
                 (ht_roof(l)/(ht_roof(l)/canyon_hwr(l)) - 0.5_r8)) * &
                 exp(-0.5_r8*canyon_hwr(l)*(1._r8-(wind_hgt_canyon(l)/ht_roof(l))))
         else  ! skimming flow
            canyon_u_wind(l) = canyontop_wind(l) * (2._r8/rpi) * &
                 exp(-0.5_r8*canyon_hwr(l)*(1._r8-(wind_hgt_canyon(l)/ht_roof(l))))
         end if

      end do

      ! Compute fluxes - Follows CLM approach for bare soils (Oleson et al 2004)

      do fl = 1, num_urbanl
         l = filter_urbanl(fl)
         g = lun%gridcell(l)

         thm_g(l) = forc_t(g) + lapse_rate*forc_hgt_t_patch(lun%patchi(l))
         thv_g(l) = forc_th(g)*(1._r8+0.61_r8*forc_q(g))
         dth(l)   = thm_g(l)-taf(l)
         dqh(l)   = forc_q(g)-qaf(l)
         dthv     = dth(l)*(1._r8+0.61_r8*forc_q(g))+0.61_r8*forc_th(g)*dqh(l)
         zldis(l) = forc_hgt_u_patch(lun%patchi(l)) - z_d_town(l)

         ! Initialize Monin-Obukhov length and wind speed including convective velocity

         call frictionvel_inst%MoninObukIni(ur(l), thv_g(l), dthv, zldis(l), z_0_town(l), um(l), obu(l))

      end do

      ! Initialize conductances
      wtus_roof(begl:endl)        = 0._r8
      wtus_road_perv(begl:endl)   = 0._r8
      wtus_road_imperv(begl:endl) = 0._r8
      wtus_sunwall(begl:endl)     = 0._r8
      wtus_shadewall(begl:endl)   = 0._r8
      wtuq_roof(begl:endl)        = 0._r8
      wtuq_road_perv(begl:endl)   = 0._r8
      wtuq_road_imperv(begl:endl) = 0._r8
      wtuq_sunwall(begl:endl)     = 0._r8
      wtuq_shadewall(begl:endl)   = 0._r8
      wtus_roof_unscl(begl:endl)        = 0._r8
      wtus_road_perv_unscl(begl:endl)   = 0._r8
      wtus_road_imperv_unscl(begl:endl) = 0._r8
      wtus_sunwall_unscl(begl:endl)     = 0._r8
      wtus_shadewall_unscl(begl:endl)   = 0._r8
      wtuq_roof_unscl(begl:endl)        = 0._r8
      wtuq_road_perv_unscl(begl:endl)   = 0._r8
      wtuq_road_imperv_unscl(begl:endl) = 0._r8
      wtuq_sunwall_unscl(begl:endl)     = 0._r8
      wtuq_shadewall_unscl(begl:endl)   = 0._r8

      ! Start stability iteration

      do iter = 1,niters

         ! Get friction velocity, relation for potential
         ! temperature and humidity profiles of surface boundary layer.

         if (num_urbanl > 0) then
            call frictionvel_inst%FrictionVelocity(begl, endl, &
                 num_urbanl, filter_urbanl, &
                 z_d_town(begl:endl), z_0_town(begl:endl), z_0_town(begl:endl), z_0_town(begl:endl), &
                 obu(begl:endl), iter, ur(begl:endl), um(begl:endl), ustar(begl:endl), &
                 temp1(begl:endl), temp2(begl:endl), temp12m(begl:endl), temp22m(begl:endl), fm(begl:endl), &
                 landunit_index=.true.)
         end if

         do fl = 1, num_urbanl
            l = filter_urbanl(fl)
            g = lun%gridcell(l)

            ! Determine aerodynamic resistance to fluxes from urban canopy air to
            ! atmosphere

            ramu(l) = 1._r8/(ustar(l)*ustar(l)/um(l))
            rahu(l) = 1._r8/(temp1(l)*ustar(l))
            rawu(l) = 1._r8/(temp2(l)*ustar(l))

            ! Determine magnitude of canyon wind by using horizontal wind determined
            ! previously and vertical wind from friction velocity (Masson 2000)

            canyon_wind(l) = sqrt(canyon_u_wind(l)**2._r8 + ustar(l)**2._r8)

            ! Determine canyon_resistance (currently this single resistance determines the
            ! resistance from urban surfaces (roof, pervious and impervious road, sunlit and
            ! shaded walls) to urban canopy air, since it is only dependent on wind speed
            ! Also from Masson 2000.

            canyon_resistance(l) = cpair * forc_rho(g) / (11.8_r8 + 4.2_r8*canyon_wind(l))

         end do

         ! This is the first term in the equation solutions for urban canopy air temperature
         ! and specific humidity (numerator) and is a landunit quantity
         do fl = 1, num_urbanl
            l = filter_urbanl(fl)
            g = lun%gridcell(l)

            taf_numer(l) = thm_g(l)/rahu(l)
            taf_denom(l) = 1._r8/rahu(l)
            qaf_numer(l) = forc_q(g)/rawu(l)
            qaf_denom(l) = 1._r8/rawu(l)

            ! First term needed for derivative of heat fluxes
            wtas(l) = 1._r8/rahu(l)
            wtaq(l) = 1._r8/rawu(l)

         end do


         ! Gather other terms for other urban columns for numerator and denominator of
         ! equations for urban canopy air temperature and specific humidity

         do fc = 1,num_urbanc
            c = filter_urbanc(fc)
            l = col%landunit(c)

            if (ctype(c) == icol_roof) then

               ! scaled sensible heat conductance
               wtus(c) = wtlunit_roof(l)/canyon_resistance(l)
               wtus_roof(l) = wtus(c)
               ! unscaled sensible heat conductance
               wtus_roof_unscl(l) = 1._r8/canyon_resistance(l)

               if (snow_depth(c) > 0._r8) then
                  fwet_roof = min(snow_depth(c)/0.05_r8, 1._r8)
               else
                  fwet_roof = (max(0._r8, h2osoi_liq(c,1)+h2osoi_ice(c,1))/pondmx_urban)**0.666666666666_r8
                  fwet_roof = min(fwet_roof,1._r8)
               end if
               if (qaf(l) > qg(c)) then 
                  fwet_roof = 1._r8
               end if
               ! scaled latent heat conductance
               wtuq(c) = fwet_roof*(wtlunit_roof(l)/canyon_resistance(l))
               wtuq_roof(l) = wtuq(c)
               ! unscaled latent heat conductance
               wtuq_roof_unscl(l) = fwet_roof*(1._r8/canyon_resistance(l))
               if ( IsSimpleBuildTemp() ) call simple_wasteheatfromac( &
               eflx_urban_ac_col(c), eflx_urban_heat_col(c), eflx_wasteheat_roof(l), &
               eflx_heat_from_ac_roof(l) )

            else if (ctype(c) == icol_road_perv) then

               ! scaled sensible heat conductance
               wtus(c) = wtroad_perv(l)*(1._r8-wtlunit_roof(l))/canyon_resistance(l)
               wtus_road_perv(l) = wtus(c)
               ! unscaled sensible heat conductance
               wtus_road_perv_unscl(l) = 1._r8/canyon_resistance(l)

               ! scaled latent heat conductance
               wtuq(c) = wtroad_perv(l)*(1._r8-wtlunit_roof(l))/canyon_resistance(l)
               wtuq_road_perv(l) = wtuq(c)
               ! unscaled latent heat conductance
               wtuq_road_perv_unscl(l) = 1._r8/canyon_resistance(l)

            else if (ctype(c) == icol_road_imperv) then

               ! scaled sensible heat conductance
               wtus(c) = (1._r8-wtroad_perv(l))*(1._r8-wtlunit_roof(l))/canyon_resistance(l)
               wtus_road_imperv(l) = wtus(c)
               ! unscaled sensible heat conductance
               wtus_road_imperv_unscl(l) = 1._r8/canyon_resistance(l)

               if (snow_depth(c) > 0._r8) then
                  fwet_road_imperv = min(snow_depth(c)/0.05_r8, 1._r8)
               else
                  fwet_road_imperv = (max(0._r8, h2osoi_liq(c,1)+h2osoi_ice(c,1))/pondmx_urban)**0.666666666666_r8
                  fwet_road_imperv = min(fwet_road_imperv,1._r8)
               end if
               if (qaf(l) > qg(c)) then 
                  fwet_road_imperv = 1._r8
               end if
               ! scaled latent heat conductance
               wtuq(c) = fwet_road_imperv*(1._r8-wtroad_perv(l))*(1._r8-wtlunit_roof(l))/canyon_resistance(l)
               wtuq_road_imperv(l) = wtuq(c)
               ! unscaled latent heat conductance
               wtuq_road_imperv_unscl(l) = fwet_road_imperv*(1._r8/canyon_resistance(l))

            else if (ctype(c) == icol_sunwall) then

               ! scaled sensible heat conductance
               wtus(c) = canyon_hwr(l)*(1._r8-wtlunit_roof(l))/canyon_resistance(l)
               wtus_sunwall(l) = wtus(c)
               ! unscaled sensible heat conductance
               wtus_sunwall_unscl(l) = 1._r8/canyon_resistance(l)

               ! scaled latent heat conductance
               wtuq(c) = 0._r8
               wtuq_sunwall(l) = wtuq(c)
               ! unscaled latent heat conductance
               wtuq_sunwall_unscl(l) = 0._r8
               if ( IsSimpleBuildTemp() ) call simple_wasteheatfromac( eflx_urban_ac_col(c),     &
                                                     eflx_urban_heat_col(c), eflx_wasteheat_sunwall(l), &
                                                     eflx_heat_from_ac_sunwall(l) )

            else if (ctype(c) == icol_shadewall) then

               ! scaled sensible heat conductance
               wtus(c) = canyon_hwr(l)*(1._r8-wtlunit_roof(l))/canyon_resistance(l)
               wtus_shadewall(l) = wtus(c)
               ! unscaled sensible heat conductance
               wtus_shadewall_unscl(l) = 1._r8/canyon_resistance(l)

               ! scaled latent heat conductance
               wtuq(c) = 0._r8
               wtuq_shadewall(l) = wtuq(c)
               ! unscaled latent heat conductance
               wtuq_shadewall_unscl(l) = 0._r8
               if ( IsSimpleBuildTemp() ) call simple_wasteheatfromac( eflx_urban_ac_col(c),     &
                                                     eflx_urban_heat_col(c), eflx_wasteheat_shadewall(l), &
                                                     eflx_heat_from_ac_shadewall(l) )

            else
               write(iulog,*) 'c, ctype, pi = ', c, ctype(c), pi
               write(iulog,*) 'Column indices for: shadewall, sunwall, road_imperv, road_perv, roof: '
               write(iulog,*) icol_shadewall, icol_sunwall, icol_road_imperv, icol_road_perv, icol_roof
               call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, &
                    msg="ERROR, ctype out of range"//errmsg(sourcefile, __LINE__))
            end if

            taf_numer(l) = taf_numer(l) + t_grnd(c)*wtus(c)
            taf_denom(l) = taf_denom(l) + wtus(c)
            qaf_numer(l) = qaf_numer(l) + qg(c)*wtuq(c)
            qaf_denom(l) = qaf_denom(l) + wtuq(c)

         end do

         ! Calculate new urban canopy air temperature and specific humidity

         call wasteheat( bounds, num_urbanl, filter_urbanl, eflx_wasteheat_roof, eflx_wasteheat_sunwall, &
                         eflx_wasteheat_shadewall, eflx_heat_from_ac_roof, eflx_heat_from_ac_sunwall,    &
                         eflx_heat_from_ac_shadewall, energyflux_inst )

         do fl = 1, num_urbanl
            l = filter_urbanl(fl)
            g = lun%gridcell(l)

            ! Calculate traffic heat flux
            ! Only comes from impervious road
            eflx_traffic(l) = (1._r8-wtlunit_roof(l))*(1._r8-wtroad_perv(l))* &
                 eflx_traffic_factor(l)

            taf(l) = taf_numer(l)/taf_denom(l)
            qaf(l) = qaf_numer(l)/qaf_denom(l)

            wts_sum(l) = wtas(l) + wtus_roof(l) + wtus_road_perv(l) + &
                 wtus_road_imperv(l) + wtus_sunwall(l) + wtus_shadewall(l)

            wtq_sum(l) = wtaq(l) + wtuq_roof(l) + wtuq_road_perv(l) + &
                 wtuq_road_imperv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)

         end do

         ! This section of code is not required if niters = 1
         ! Determine stability using new taf and qaf
         ! TODO: Some of these constants replicate what is in FrictionVelocity and BareGround fluxes should consildate. EBK
         do fl = 1, num_urbanl
            l = filter_urbanl(fl)
            g = lun%gridcell(l)

            dth(l) = thm_g(l)-taf(l)
            dqh(l) = forc_q(g)-qaf(l)
            tstar = temp1(l)*dth(l)
            qstar = temp2(l)*dqh(l)
            thvstar = tstar*(1._r8+0.61_r8*forc_q(g)) + 0.61_r8*forc_th(g)*qstar
            zeta = zldis(l)*vkc*grav*thvstar/(ustar(l)**2*thv_g(l))

            if (zeta >= 0._r8) then                   !stable
               zeta = min(zetamax,max(zeta,0.01_r8))
               um(l) = max(ur(l),0.1_r8)
            else                                      !unstable
               zeta = max(-100._r8,min(zeta,-0.01_r8))
               wc = beta(l)*(-grav*ustar(l)*thvstar*zii(l)/thv_g(l))**0.333_r8
               um(l) = sqrt(ur(l)*ur(l) + wc*wc)
            end if

            obu(l) = zldis(l)/zeta
         end do

      end do   ! end iteration

      ! Determine fluxes from canyon surfaces

      ! the following initializations are needed to ensure that the values are 0 over non-
      ! active urban Patches
      eflx_sh_grnd_scale(bounds%begp : bounds%endp) = 0._r8
      qflx_evap_soi_scale(bounds%begp : bounds%endp) = 0._r8

      do f = 1, num_urbanp

         p = filter_urbanp(f)
         c = patch%column(p)
         g = patch%gridcell(p)
         l = patch%landunit(p)

         ram1(p) = ramu(l)  !pass value to global variable

         ! Upward and downward canopy longwave are zero

         ulrad(p)  = 0._r8
         dlrad(p)  = 0._r8

         ! Derivative of sensible and latent heat fluxes with respect to 
         ! ground temperature

         if (ctype(c) == icol_roof) then
            cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_road_perv(l) +  &
                 wtus_road_imperv(l) + wtus_sunwall(l) + wtus_shadewall(l)) * &
                 (wtus_roof_unscl(l)/wts_sum(l))
            cgrndl(p) = forc_rho(g) * (wtaq(l) + wtuq_road_perv(l) +  &
                 wtuq_road_imperv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)) * &
                 (wtuq_roof_unscl(l)/wtq_sum(l))*dqgdT(c)
         else if (ctype(c) == icol_road_perv) then
            cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) +  &
                 wtus_road_imperv(l) + wtus_sunwall(l) + wtus_shadewall(l)) * &
                 (wtus_road_perv_unscl(l)/wts_sum(l))
            cgrndl(p) = forc_rho(g) * (wtaq(l) + wtuq_roof(l) +  &
                 wtuq_road_imperv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)) * &
                 (wtuq_road_perv_unscl(l)/wtq_sum(l))*dqgdT(c)
         else if (ctype(c) == icol_road_imperv) then
            cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) +  &
                 wtus_road_perv(l) + wtus_sunwall(l) + wtus_shadewall(l)) * &
                 (wtus_road_imperv_unscl(l)/wts_sum(l))
            cgrndl(p) = forc_rho(g) * (wtaq(l) + wtuq_roof(l) +  &
                 wtuq_road_perv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)) * &
                 (wtuq_road_imperv_unscl(l)/wtq_sum(l))*dqgdT(c)
         else if (ctype(c) == icol_sunwall) then
            cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) +  &
                 wtus_road_perv(l) + wtus_road_imperv(l) + wtus_shadewall(l)) * &
                 (wtus_sunwall_unscl(l)/wts_sum(l))
            cgrndl(p) = 0._r8
         else if (ctype(c) == icol_shadewall) then
            cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) +  &
                 wtus_road_perv(l) + wtus_road_imperv(l) + wtus_sunwall(l)) * &
                 (wtus_shadewall_unscl(l)/wts_sum(l))
            cgrndl(p) = 0._r8
         end if
         cgrnd(p)  = cgrnds(p) + cgrndl(p)*htvp(c)

         ! Surface fluxes of momentum, sensible and latent heat

         taux(p)          = -forc_rho(g)*forc_u(g)/ramu(l)
         tauy(p)          = -forc_rho(g)*forc_v(g)/ramu(l)

         ! Use new canopy air temperature
         dth(l) = taf(l) - t_grnd(c)

         if (ctype(c) == icol_roof) then
            eflx_sh_grnd(p)  = -forc_rho(g)*cpair*wtus_roof_unscl(l)*dth(l)
            eflx_sh_snow(p)  = 0._r8
            eflx_sh_soil(p)  = 0._r8
            eflx_sh_h2osfc(p)= 0._r8
         else if (ctype(c) == icol_road_perv) then
            eflx_sh_grnd(p)  = -forc_rho(g)*cpair*wtus_road_perv_unscl(l)*dth(l)
            eflx_sh_snow(p)  = 0._r8
            eflx_sh_soil(p)  = 0._r8
            eflx_sh_h2osfc(p)= 0._r8
         else if (ctype(c) == icol_road_imperv) then
            eflx_sh_grnd(p)  = -forc_rho(g)*cpair*wtus_road_imperv_unscl(l)*dth(l)
            eflx_sh_snow(p)  = 0._r8
            eflx_sh_soil(p)  = 0._r8
            eflx_sh_h2osfc(p)= 0._r8
         else if (ctype(c) == icol_sunwall) then
            eflx_sh_grnd(p)  = -forc_rho(g)*cpair*wtus_sunwall_unscl(l)*dth(l)
            eflx_sh_snow(p)  = 0._r8
            eflx_sh_soil(p)  = 0._r8
            eflx_sh_h2osfc(p)= 0._r8
         else if (ctype(c) == icol_shadewall) then
            eflx_sh_grnd(p)  = -forc_rho(g)*cpair*wtus_shadewall_unscl(l)*dth(l)
            eflx_sh_snow(p)  = 0._r8
            eflx_sh_soil(p)  = 0._r8
            eflx_sh_h2osfc(p)= 0._r8
         end if

         eflx_sh_tot(p)   = eflx_sh_grnd(p)
         eflx_sh_tot_u(p) = eflx_sh_tot(p)

         dqh(l) = qaf(l) - qg(c)

         if (ctype(c) == icol_roof) then
            qflx_evap_soi(p) = -forc_rho(g)*wtuq_roof_unscl(l)*dqh(l)
         else if (ctype(c) == icol_road_perv) then
            ! Evaporation assigned to soil term if dew or snow
            ! or if no liquid water available in soil column
            if (dqh(l) > 0._r8 .or. frac_sno(c) > 0._r8 .or. soilalpha_u(c) <= 0._r8) then
               qflx_evap_soi(p) = -forc_rho(g)*wtuq_road_perv_unscl(l)*dqh(l)
               qflx_tran_veg(p) = 0._r8
               ! Otherwise, evaporation assigned to transpiration term
            else
               qflx_evap_soi(p) = 0._r8
               qflx_tran_veg(p) = -forc_rho(g)*wtuq_road_perv_unscl(l)*dqh(l)
            end if
            qflx_evap_veg(p) = qflx_tran_veg(p)
         else if (ctype(c) == icol_road_imperv) then
            qflx_evap_soi(p) = -forc_rho(g)*wtuq_road_imperv_unscl(l)*dqh(l)
         else if (ctype(c) == icol_sunwall) then
            qflx_evap_soi(p) = 0._r8
         else if (ctype(c) == icol_shadewall) then
            qflx_evap_soi(p) = 0._r8
         end if

         ! SCALED sensible and latent heat flux for error check
         eflx_sh_grnd_scale(p)  = -forc_rho(g)*cpair*wtus(c)*dth(l)
         qflx_evap_soi_scale(p) = -forc_rho(g)*wtuq(c)*dqh(l)

      end do

      ! Check to see that total sensible and latent heat equal the sum of
      ! the scaled heat fluxes above
      do fl = 1, num_urbanl
         l = filter_urbanl(fl)
         g = lun%gridcell(l)
         eflx(l)       = -(forc_rho(g)*cpair/rahu(l))*(thm_g(l) - taf(l))
         qflx(l)       = -(forc_rho(g)/rawu(l))*(forc_q(g) - qaf(l))
         eflx_scale(l) = sum(eflx_sh_grnd_scale(lun%patchi(l):lun%patchf(l)))
         qflx_scale(l) = sum(qflx_evap_soi_scale(lun%patchi(l):lun%patchf(l)))
         eflx_err(l)   = eflx_scale(l) - eflx(l)
         qflx_err(l)   = qflx_scale(l) - qflx(l)
      end do

      found = .false.
      do fl = 1, num_urbanl
         l = filter_urbanl(fl)
         if (abs(eflx_err(l)) > 0.01_r8) then
            found = .true.
            indexl = l
            exit
         end if
      end do
      if ( found ) then
         write(iulog,*)'WARNING:  Total sensible heat does not equal sum of scaled heat fluxes for urban columns ',&
              ' nstep = ',nstep,' indexl= ',indexl,' eflx_err= ',eflx_err(indexl)
         if (abs(eflx_err(indexl)) > .01_r8) then
            write(iulog,*)'clm model is stopping - error is greater than .01 W/m**2'
            write(iulog,*)'eflx_scale    = ',eflx_scale(indexl)
            write(iulog,*)'eflx_sh_grnd_scale: ',eflx_sh_grnd_scale(lun%patchi(indexl):lun%patchf(indexl))
            write(iulog,*)'eflx          = ',eflx(indexl)
            call endrun(subgrid_index=indexl, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__))
         end if
      end if

      found = .false.
      do fl = 1, num_urbanl
         l = filter_urbanl(fl)
         ! 4.e-9 kg/m**2/s = 0.01 W/m**2
         if (abs(qflx_err(l)) > 4.e-9_r8) then
            found = .true.
            indexl = l
            exit
         end if
      end do
      if ( found ) then
         write(iulog,*)'WARNING:  Total water vapor flux does not equal sum of scaled water vapor fluxes for urban columns ',&
              ' nstep = ',nstep,' indexl= ',indexl,' qflx_err= ',qflx_err(indexl)
         if (abs(qflx_err(indexl)) > 4.e-9_r8) then
            write(iulog,*)'clm model is stopping - error is greater than 4.e-9 kg/m**2/s'
            write(iulog,*)'qflx_scale    = ',qflx_scale(indexl)
            write(iulog,*)'qflx          = ',qflx(indexl)
            call endrun(subgrid_index=indexl, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__))
         end if
      end if

      ! Gather terms required to determine internal building temperature

      if ( IsSimpleBuildTemp() ) call calc_simple_internal_building_temp( &
                  bounds, num_urbanc, filter_urbanc, num_urbanl, filter_urbanl,     &
                  temperature_inst)

      ! No roots for urban except for pervious road

      do j = 1, nlevgrnd
         do f = 1, num_urbanp
            p = filter_urbanp(f)
            c = patch%column(p)
            if (ctype(c) == icol_road_perv) then
               rootr(p,j) = rootr_road_perv(c,j)
            else
               rootr(p,j) = 0._r8
            end if
         end do
      end do

      do f = 1, num_urbanp

         p = filter_urbanp(f)
         c = patch%column(p)
         g = patch%gridcell(p)
         l = patch%landunit(p)

         ! Use urban canopy air temperature and specific humidity to represent 
         ! 2-m temperature and humidity

         t_ref2m(p) = taf(l)
         q_ref2m(p) = qaf(l)
         t_ref2m_u(p) = taf(l)

         ! 2 m height relative humidity

         call QSat(t_ref2m(p), forc_pbot(g), qsat_ref2m, &
              es = e_ref2m)
         rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8)
         rh_ref2m_u(p) = rh_ref2m(p)

         ! Human Heat Stress
         if ( all_human_stress_indices .or. fast_human_stress_indices )then
            call KtoC(t_ref2m(p), tc_ref2m(p))
            call VaporPres(rh_ref2m(p), e_ref2m, vap_ref2m(p))
            call Wet_BulbS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p))
            call HeatIndex(tc_ref2m(p), rh_ref2m(p), nws_hi_ref2m(p))
            call AppTemp(tc_ref2m(p), vap_ref2m(p), u10_clm(p), appar_temp_ref2m(p))
            call swbgt(tc_ref2m(p), vap_ref2m(p), swbgt_ref2m(p))
            call hmdex(tc_ref2m(p), vap_ref2m(p), humidex_ref2m(p))
            call dis_coiS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p), discomf_index_ref2mS(p))
            if ( all_human_stress_indices ) then
               call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(g), rh_ref2m(p), q_ref2m(p), &
                             teq_ref2m(p), ept_ref2m(p), wb_ref2m(p))
               call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p))
               call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p))
               call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p))
            end if
  
            wbt_ref2m_u(p)            = wbt_ref2m(p)
            nws_hi_ref2m_u(p)         = nws_hi_ref2m(p)
            appar_temp_ref2m_u(p)     = appar_temp_ref2m(p)
            swbgt_ref2m_u(p)          = swbgt_ref2m(p)
            humidex_ref2m_u(p)        = humidex_ref2m(p)
            discomf_index_ref2mS_u(p) = discomf_index_ref2mS(p)
            if ( all_human_stress_indices ) then
               teq_ref2m_u(p)            = teq_ref2m(p)
               ept_ref2m_u(p)            = ept_ref2m(p)
               wb_ref2m_u(p)             = wb_ref2m(p)
               discomf_index_ref2m_u(p)  = discomf_index_ref2m(p)
               thic_ref2m_u(p)           = thic_ref2m(p)
               thip_ref2m_u(p)           = thip_ref2m(p)
               swmp80_ref2m_u(p)         = swmp80_ref2m(p)
               swmp65_ref2m_u(p)         = swmp65_ref2m(p)
            end if
         end if

         ! Variables needed by history tape

         t_veg(p) = forc_t(g)

      end do

    end associate

  end subroutine UrbanFluxes

  !----------------------------------------------------------------------- 
  !BOP
  !
  ! !IROUTINE: wasteheat
  !
  ! !INTERFACE:
  subroutine wasteheat( bounds, num_urbanl, filter_urbanl, eflx_wasteheat_roof, eflx_wasteheat_sunwall, &
                        eflx_wasteheat_shadewall, eflx_heat_from_ac_roof, eflx_heat_from_ac_sunwall,    &
                        eflx_heat_from_ac_shadewall, energyflux_inst )
    ! !DESCRIPTION:
    !
    ! Calculate the wasteheat flux from urban heating or air-conditioning.
    !
    ! !USES:
    use clm_varcon         , only : ht_wasteheat_factor, ac_wasteheat_factor, &
                                    wasteheat_limit
    use EnergyFluxType     , only : energyflux_type
    use UrbanParamsType    , only : IsProgBuildTemp
    implicit none
    ! !ARGUMENTS:
    type(bounds_type), intent(in) :: bounds  ! bounds
    integer , intent(in) :: num_urbanl       ! number of urban landunits in clump
    integer , intent(in) :: filter_urbanl(:) ! urban landunit filter
    real(r8)            , intent(in)  :: eflx_wasteheat_roof(bounds%begl:bounds%endl)
    real(r8)            , intent(in)  :: eflx_wasteheat_sunwall(bounds%begl:bounds%endl)
    real(r8)            , intent(in)  :: eflx_wasteheat_shadewall(bounds%begl:bounds%endl)
    real(r8)            , intent(in)  :: eflx_heat_from_ac_roof(bounds%begl:bounds%endl)
    real(r8)            , intent(in)  :: eflx_heat_from_ac_sunwall(bounds%begl:bounds%endl)
    real(r8)            , intent(in)  :: eflx_heat_from_ac_shadewall(bounds%begl:bounds%endl)
    type(energyflux_type) , intent(inout)  :: energyflux_inst  ! data on landunit energy flux

    ! !LOCAL VARIABLES:
    integer fl, l, g
    !EOP
    !----------------------------------------------------------------------- 

    associate(&
     lgridcell        => lun%gridcell     , & ! Input:  [integer (:)    ]  gridcell of corresponding landunit                 
     canyon_hwr       => lun%canyon_hwr   , & ! Input:  [real(r8) (:)]    ratio of building height to street width 
     wtlunit_roof     => lun%wtlunit_roof , & ! Input:  [real(r8) (:)]    weight of roof with respect to landunit
     eflx_wasteheat   => energyflux_inst%eflx_wasteheat_lun    , & ! Output:  [real(r8) (:)]  sensible heat flux from urban heating/cooling sources of waste heat (W/m**2)
     eflx_heat_from_ac=> energyflux_inst%eflx_heat_from_ac_lun , & ! Output:  [real(r8) (:)]  sensible heat flux put back into canyon due to removal by AC (W/m**2)
     eflx_urban_ac    => energyflux_inst%eflx_urban_ac_lun     , & ! Input:  [real(r8) (:)]  urban air conditioning flux (W/m**2)              
     eflx_urban_heat  => energyflux_inst%eflx_urban_heat_lun     & ! Input:  [real(r8) (:)]  urban heating flux (W/m**2)                       
    )
    do fl = 1, num_urbanl
       l = filter_urbanl(fl)
       g = lgridcell(l)
       if      ( IsSimpleBuildTemp() )then
          ! Total waste heat and heat from AC is sum of heat for walls and roofs
          ! accounting for different surface areas
          eflx_wasteheat(l) = wtlunit_roof(l)*eflx_wasteheat_roof(l) + &
             (1._r8-wtlunit_roof(l))*(canyon_hwr(l)*(eflx_wasteheat_sunwall(l) + &
             eflx_wasteheat_shadewall(l)))

       else if ( IsProgBuildTemp() )then
          ! wasteheat from heating/cooling
          if (trim(urban_hac) == urban_wasteheat_on) then
            eflx_wasteheat(l) = ac_wasteheat_factor * eflx_urban_ac(l) + &
                                ht_wasteheat_factor * eflx_urban_heat(l)
          else
            eflx_wasteheat(l) = 0._r8
          end if
       end if

       ! Limit wasteheat to ensure that we don't get any unrealistically strong 
       ! positive feedbacks due to AC in a warmer climate
       eflx_wasteheat(l) = min(eflx_wasteheat(l),wasteheat_limit)

       if      ( IsSimpleBuildTemp() )then
          eflx_heat_from_ac(l) = wtlunit_roof(l)*eflx_heat_from_ac_roof(l) + &
          (1._r8-wtlunit_roof(l))*(canyon_hwr(l)*(eflx_heat_from_ac_sunwall(l) + &
          eflx_heat_from_ac_shadewall(l)))

       else if ( IsProgBuildTemp() )then
          ! If air conditioning on, always replace heat removed with heat into canyon
          if (trim(urban_hac) == urban_hac_on .or. trim(urban_hac) == urban_wasteheat_on) then
            eflx_heat_from_ac(l) = abs(eflx_urban_ac(l))
          else
            eflx_heat_from_ac(l) = 0._r8
          end if
       end if
    end do
  end associate
  end subroutine wasteheat

  !----------------------------------------------------------------------- 
  !BOP
  !
  ! !IROUTINE: simple_wasteheatfromac
  !
  ! !INTERFACE:
  subroutine simple_wasteheatfromac( eflx_urban_ac, eflx_urban_heat, eflx_wasteheat, &
                  eflx_heat_from_ac )
    !----------------------------------------------------------------------- 
    ! !DESCRIPTION: 
    !
    ! Calculate waste heat from Air conditioning with the simpler method introduced
    ! in CLM4.5.
    !
    ! !USES:
    use clm_varcon         , only : ht_wasteheat_factor, ac_wasteheat_factor
    implicit none
    ! !ARGUMENTS:
    real(r8), intent(in)  :: eflx_urban_ac
    real(r8), intent(in)  :: eflx_urban_heat
    real(r8), intent(out) :: eflx_wasteheat
    real(r8), intent(out) :: eflx_heat_from_ac

    ! wasteheat from heating/cooling
    if (trim(urban_hac) == urban_wasteheat_on) then
       eflx_wasteheat = ac_wasteheat_factor * eflx_urban_ac + &
            ht_wasteheat_factor * eflx_urban_heat
    else
       eflx_wasteheat = 0._r8
    end if

    ! If air conditioning on, always replace heat removed with heat into canyon
    if (trim(urban_hac) == urban_hac_on .or. trim(urban_hac) == urban_wasteheat_on) then
       eflx_heat_from_ac = abs(eflx_urban_ac)
    else
       eflx_heat_from_ac = 0._r8
    end if

  end subroutine simple_wasteheatfromac

  !----------------------------------------------------------------------- 
  !BOP
  !
  ! !IROUTINE: calc_simple_internal_building_temp
  !
  ! !INTERFACE:
  subroutine calc_simple_internal_building_temp( bounds, num_urbanc, filter_urbanc, &
                  num_urbanl, filter_urbanl, temperature_inst )
  !----------------------------------------------------------------------- 
  ! !DESCRIPTION: 
  !
  ! Calculate the internal building temperature, based on the simpler method introduced
  ! in CLM4.5.
  !
  ! !USES:
    use clm_varpar     , only : nlevurb
    use column_varcon  , only : icol_roof, icol_sunwall, icol_shadewall
    use LandunitType   , only : landunit_type
    use ColumnType     , only : column_type
    use TemperatureType, only : temperature_type

    implicit none
  ! !ARGUMENTS:
    type(bounds_type), intent(in) :: bounds           ! bounds
    integer          , intent(in) :: num_urbanl       ! number of urban landunits in clump
    integer          , intent(in) :: filter_urbanl(:) ! urban landunit filter
    integer          , intent(in) :: num_urbanc       ! number of urban columns in clump
    integer          , intent(in) :: filter_urbanc(:) ! urban column filter
    type(temperature_type), intent(inout)  :: temperature_inst ! temperature variables
  ! !LOCAL VARIABLES:
    ! Gather terms required to determine internal building temperature
    integer  :: fl,fc,l,c                                    ! indices
    real(r8) :: t_sunwall_innerl(bounds%begl:bounds%endl)    ! temp of inner layer of sunwall (K)
    real(r8) :: t_shadewall_innerl(bounds%begl:bounds%endl)  ! temp of inner layer of shadewall (K)
    real(r8) :: t_roof_innerl(bounds%begl:bounds%endl)       ! temp of inner layer of roof (K)
    real(r8) :: lngth_roof                                   ! length of roof (m)
  !EOP
  !----------------------------------------------------------------------- 

    associate(&
     t_soisno      =>    temperature_inst%t_soisno_col  , & ! Input:  [real(r8) (:,:)]  soil temperature (K) 
     ht_roof       =>    lun%ht_roof                    , & ! Input:  [real(r8) (:)]    height of urban roof (m)
     canyon_hwr    =>    lun%canyon_hwr                 , & ! Input:  [real(r8) (:)]    ratio of building height to street width 
     wtlunit_roof  =>    lun%wtlunit_roof               , & ! Input:  [real(r8) (:)]    weight of roof with respect to landunit
     t_building    =>    temperature_inst%t_building_lun  & ! Output: [real(r8) (:)]  internal building temperature (K) 
    )

    do fc = 1,num_urbanc
       c = filter_urbanc(fc)
       l = col%landunit(c)

       if      (col%itype(c) == icol_roof     ) then
          t_roof_innerl(l)      = t_soisno(c,nlevurb)
       else if (col%itype(c) == icol_sunwall  ) then
          t_sunwall_innerl(l)   = t_soisno(c,nlevurb)
       else if (col%itype(c) == icol_shadewall) then
          t_shadewall_innerl(l) = t_soisno(c,nlevurb)
       end if

    end do

    ! Calculate internal building temperature
    do fl = 1, num_urbanl
       l = filter_urbanl(fl)
     
       lngth_roof = (ht_roof(l)/canyon_hwr(l))*wtlunit_roof(l)/(1._r8-wtlunit_roof(l))
       t_building(l) = (ht_roof(l)*(t_shadewall_innerl(l) + t_sunwall_innerl(l)) &
                       +lngth_roof*t_roof_innerl(l))/(2._r8*ht_roof(l)+lngth_roof)
    end do

  end associate

  end subroutine calc_simple_internal_building_temp

  !----------------------------------------------------------------------- 

end module UrbanFluxesMod
